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Abstract 

We study a Josephson junction ladder in a magnetic field in the absence 

of charging effects via a transfer matrix formalism. The eigenvalues of the 

transfer matrix are found numerically, giving a determination of the different 

phases of the ladder. The spatial periodicity of the ground state exhibits 

a devil's staircase as a function of the magnetic flux filling factor /. If the 

transverse Josephson coupling is varied a continuous superconducting-normal 

transition in the transverse direction is observed, analogous to the breakdown 

of the KAM trajectories in dynamical systems. We also examine how these 

properties may be affected by a current injected along the ladder. 
74.50.+r, 05.20.-y, 64.70.Rh 
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I. INTRODUCTION 



Two-dimensional arrays of Josephson junctions have attracted much recent theoretical 
and experimental attention M. Interesting physics arises as a result of competing vortex- 
vortex and vortex-lattice interactions. It is also considered to be a convenient experimental 
realization of the frustrated XY models. In this paper, we expand and elaborate on our 
previous letter |2| on the simplest of such system, namely the Josephson junction ladder 
( JJL) M shown in Figure |. 

To construct the system, superconducting elements are placed at the ladder sites. Below 
the bulk superconducting-normal transition temperature, the state of each element is de- 
scribed by its charge and the phase of the superconducting wave function || . In this paper we 
neglect charging effects, which corresponds to the condition that 4e 2 /C <C J, with C being 
the capacitance of the element and J the Josephson coupling. Let 9j (#'•) denote the phase on 
the upper (lower) branch of the ladder at the j'th rung. The Hamiltonian for the array J?| can 
be written in terms the gauge invariant phase differences, jj = 9j — Oj-i — (27r/0 o ) JJ_i A x dx, 
ii = Q 'i ~ °U ~ ( 27T /^o) If -i A x dx, and a, = Q) - 9 3 - (2tt/0 o ) jf A y dx: 

Ti = — y^X^x cos 7j + Jx cos Jj + J y COS (Xj), (1) 
j 

where A x and A y are the components of the magnetic vector potential along and transverse 
to the ladder, respectively, and (f> the flux quantum. The sum of the phase differences 
around a plaquette is constrained by 

7j-^ + a i -a j _i = 27r(/-n i ), (2) 

where rij = 0, ±1, ±2, ■ • • is the vortex occupancy number and / = <p/4>o with being the 
magnetic flux through a plaquette. With this constraint, it is convenient to write Eq. (]T|) in 
the form 

Ti = — Jj^{2 cosrjj cos[(a,_i — <x,)/2 + ir(f — rij)] 
j 

+ J t cos atj}, (3) 
2 



where rjj = (jj + Jj)/2, J = J x and J t = J y /J x - The Hamiltonian is symmetric under 
/ — *■ f + 1 with rij — > rij + 1, and / — > —f with rij — > —Tij, thus it is sufficient to 
study only the region < / < 0.5. Since in one dimension ordered phases occur only at zero 
temperature, the main interest is in the ground states of the ladder and the low temperature 
excitations. Note that in Eq. (|3|) rjj decouples from aj and rij, so that all the ground states 
have rjj = to minimize TC. The ground states will be among the solutions to the current 
conservation equations dTC/daj = 0: 

Jtsinatj = sin[(aj_i — Oj)/2 + n(f — rij)} 

- sin[(a j - a j+1 )/2 + n(f - n j+1 )]. (4) 

For any given / there are a host of solutions to Eq. (|]). The solution that minimizes the 
energy must be selected to obtain the ground state. 

If one expands the cos[(aj_i — Qj)/2 + n(f — rij)} term in Eq. (||) about it's maximum 
and set rjj = 0, the discrete sine-Gordan model (DSG) is obtained: 

H = - J£{k«i-i - aj)/2 + vr/] 2 + J t cos a j}, (5) 

A vortex (rij = 1) in the JJL corresponds to a kink in the DSG (in the DSG the a absorb 
the rij and are no longer restricted to (— tt, 7r]). Kardar || used this analogy to argue that 
this system should show similar behavior to the DSG which has been studied by several 
authors |8HT0|1 . This analogy is only valid for J t very small so that the inter-plaquette term 



dominates the behavior of the system making the expansion about its maximum a reasonable 
assumption. However, much of the interesting behavior of the DSG occurs in regions of large 
Jt (Jt ~ 1)- Furthermore, much of the work by Aubry || on the DSG relies on the convexity 
of the coupling potential which we do not have in the JJL. 

In this paper we formulate the problem in terms of a transfer matrix obtained from 
the full partition function of the ladder. The eigenvalues and eigenfunctions of the transfer 
matrix are found numerically to determine the phases of the ladder as functions of / and 
J t - We find that the spatial periodicity of the ground states goes through a devil's staircase 



as a function of /. We then study the properties of various ground states and the low 
temperature excitations. As J t is varied, all incommensurate ground states are found to 
undergo a superconducting-normal transition at certain J t which depends on /. In the last 
section we discuss the effects of a current. 



II. TRANSFER MATRIX FORMULATION 

The partition function for the ladder, with periodic boundary conditions is 

N w 

Z = Y[ datidrji exp{K (2 cos rjiCos[(ai-i - + 7r(/ - m)] + J t cosctj)} . (6) 

where K = J/ksT. The rji can be integrated out and n» summed over, resulting in a simple 
transfer matrix formalism for the partition function involving only the transverse phase 
differences: Z = dociP{oii-i, a*) = Tr P N . The transfer matrix elements P(a,a') are 

P{a,a') =4nexp[KJ t {cosa + cosa')/2]I {2Kcos[{a-a')/2 + iTf}), (7) 

where Iq is the zeroth order modified Bessel function = - Jq exp (x cos r))dr)). Note 

that the elements of P are real and positive, so that its largest eigenvalue Ao is real, positive 
and non-degenerate. However, since P is not symmetric (except for / = and / = 1/2) 
other eigenvalues can form complex conjugate pairs. As we will see from the correlation 
function, these complex eigenvalues determine the spatial periodicity of the ground states. 
The two point correlation function of a/s is 



D i(a -c»i) 



Kao-aO) = lim 



N— >oo 

I 



= E^(r)' (8) 

where the \ n are the eigenvalues (|A n | > |A n+ i| and n = 0, 1,2, •••), and the constants 
c n = fZj; dai'ipQ {ot!)e %a ifj^ (a') daijj^(a)e~ ia ^pQ(a). (Note that since P is not symmetric 
both right ip^ and left eigenf unctions are needed.) If Ai is real and |Ai| > |A 2 |, Eq. @ 
simplifies for large I to 

<e i(«o-a,)) = Co + c ; |Ai| > |A 2 |. 
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If Ai = A2 = |Ai|e Ji7r =, Eq. d) for large / is 

A ' 

(e^ a °- a ^) = c + ( Cl e i2nSl + c 2 e" M ) -± 
v ' A 

Note that while the correlation length is given by £ = [In | Ao/Ai|] _1 the quantity 5 = 
v4rg(Ai)/27r determines the spatial periodicity of the state. For example, Figure ^| shows 
the full correlation function (not just the first two terms) for a three periodic state, H = 
1/3. We see that the correlation function has three branches corresponding to (e^" 0-03 ' 1 -*), 

( e i(a -a 3n+1 )^ and ( e i(a -a 3n+2 )} w j th n = 0, 1, 2, • • •. 

It is fairly easy to find the eigenvalues and eigenvectors of the transfer matrix numer- 
ically to very high accuracy (see the Appendix). For / smaller than a critical value / c i 
which depends on J t , we find that both Ai and A2 are real. These two eigenvalues become 
degenerate at f c i, and then bifurcate into a complex conjugate pair. H as a function of / is 
shown in Figure || for several different values of Jt. The shape of these curves is generally 
referred to as a devil's staircase. The steps of the staircase are at S = p/q, where p and q 
are integers. These are commensurate states with p vortices in each unit cell which consists 
of q plaquettes. For small Jt, the flat steps are connected by fairly smooth curves; most 
states on the H — / curve are incommensurate states. As J t increases, more and more steps 
appear and grow at the expense of the smooth regions. It appears that at Jt — J t c ~ 0.5 
the staircase becomes complete, i.e. there is a step for every rational S and the set of / 
which correspond to irrational H has zero measure. For J t > J£, the staircase becomes 
over-complete, i.e. steps of lower order rationals grow and steps of higher order rationals 
disappear |IT[] . A phase diagram can be constructed with the phase boundaries at the step 



edg shown in Figure El. 

Another important characterization of a state is the phase density p(ot): p(a)da is the 
average fraction of all sites in the ladder with a < ai < a + da. If p(a) is a smooth function 
for a G (— 7T, 7r] at T = 0, the ground state energy is invariant under an adiabatic change 
of a's. Consequently, there is no phase coherence between upper and lower branches of the 
ladder and hence no superconductivity in the transverse direction. In this case, we say that 



the a's are unpinned. If there exist finite intervals of a on which p(a) = 0, the a's are 
pinned and there will be phase coherence between the upper and lower branches. In terms 
of the transfer matrix, the phase density is the product of the left and right eigenfunctions 
of Ac ^,p(a)=^(a)^(a). 



III. COMMENSURATE STATES 



We first discuss the case where f < f c i- These are the "Meissner" states in the sense that 
there are no vortices (n« = 0) in the ladder. The ground state is simply on = 0, 7^ = 71"/ and 
7^ = —7r/, so that there is a global "screening" current ±J x sin7r/ in the upper and lower 
branches of the ladder ||. The phase density p(a) = 5(a). The properties of the Meissner 
state can be studied by expanding Eq. (|3|) around a, = 0: 



H M = (J/4) ^[cos( 7 r/)(a i _ 1 - a,) 2 + 2J t a] 



(9) 



At finite temperature, p(a) = 5(a) peaks become thermally broadened. The fluctuations 
about «j = in the J t a 2 part of the energy of a single plaquette will be of order k{T (from 
equipartition). Hence 5ai ~ y/k^T is an estimate of the p(a) peak width. The current 
conservation Eq. (JD becomes 



= 2 (1 + J t j cos 71"/) aj — 



(10) 



Besides the ground state aj = 0, there are other two linearly independent solutions aj 
e ±i/?M Q f Eq flipp which describe collective fluctuations about the ground state, where 



Jt 



COS 7T , 



2J t 



f ~*~ \ COS7T/ ^ ^COS7r/ 



Jt 



(11) 



£m is the low temperature correlation length for the Meissner state. (Note that £m < 1 for 
J t ~ 1 making a continuum approximation invalid.) A comparison of £m to the correlation 
length obtained from the transfer matrix (or from an exact solution of the Gaussian model 



) partition function |L3|]) shows that the two are indistinguishable for k^T ' / J < 0.005. Note 



that for finite Jt this model has no zero temperature phase transition in the sense that the 
correlation length remains finite as T — > 0. Furthermore, the correlation length diverges like 
1 / y/2Jt as Jt — > 0. (or more accurately, approaches the rigid rotator value of 1/kbT, as we 
shall see below.) 

As / increases, the Meissner state becomes unstable to the formation of vortices. A 
vortex is constructed by patching the two vortex free solutions of the Eq. (|10( ) together 
anti-symmetrically about the plaquette with the vortex. Using Eq. @ for the matching 
condition, 

- sin(a* + Tr/) - sin((e 1/5M - l)a./2 + nf) + sin a* = 0. 

where ±a* is the the value of aj on the plaquette enclosing the vortex. Expanded to second 
order, this gives 

-4sin(7r/) + [2 - (e 1/?JU + 1) cos(vr/)]a, + [1 + (e 1/?JU - l) 2 ] sin(7r/)a 2 = 

To first order, the energy e v of a single vortex is found to be 

e v = 4cos(7r/) — 3a*sin(7r/) 

The zero of determines f c \. This gives f c \ ~ 0.26 for J t = 1. Extending the calculation of 
e v to second order gives f& ~ 0.28 for J t = 1 which is in good agreement with the numerical 
result from the transfer matrix. 

For f > fdi £v is negative and vortices are spontaneously created. When vortices are far 
apart their interaction is caused only by the exponentially small overlap. The corresponding 
repulsion energy is of the order Jexp(— 1/£,m)j where I is the distance between vortices. At 
finite temperatures and low densities, the vortices will be able to overcome a weak pinning 
potential and will move around fairly freely in the I sites separating the vortices. This leads 
to a free energy per plaquette of 

F = e v /l + Jcexp(-l/£ M )/l - k B T\nl/l, 

where c is a constant of order unity. Minimizing this free energy as a function of / and putting 
in the dependence of e v near f cl gives (rij) = I' 1 = [6/ In |/ c i-/|] _1 , for (f-f c i)/k B T » I. 
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We now discuss the commensurate vortex states, taking the one with 5 = 1/2 as an 
example. This state has many similarities to the Meissner state but some important differ- 
ences. The ground state is 

(a ,n ) = (arctan[(2/J t )sin(vr/)] ,0), (ai,n x ) = (-a , 1); (a<±2, n i±2 ) = (a.,^); 

so that there is a global screening current in the upper and lower branches of the ladder of 



±2ttJ(/-1/2)/- v /4 + J? and the energy is yjl + (2/J t ) 2 sm 2 (nf). Global screening;, which is 
absent in an infinite 2D array, is the key reason for the existence of the steps at 5 = pj q. It is 
easy to see that the symmetry of this 5 = 1/2 vortex state is that of the (antiferromagnetic) 
Ising model. The low temperature excitations are domain boundaries between the two 
degenerate ground states. The energy of the domain boundary Je^ can be estimated using 
similar methods to those used to derive e v for the Meissner state. We found that £& = 



e° — (it 2 1 J 4 + J 2 )\f — 1/2|, where e° depends only on J t . Thus the correlation length 
diverges with temperature as £ ~ exp(2 Je^/Zc^T). The transition from the 5 = 1/2 state 
to nearby vortex states happens when / is such that et, = 0; it is similar to the transition 
from the Meissner state to its nearby vortex states. All other steps 5 = p/q can be analyzed 
similarly. For comparison, we have evaluated £ for various values of / and T from the 
transfer matrix and found that £ fits £ ~ exp(2Jeb/kBT) (typically over several decades) 
at low temperature. The value of e& as a function of / is shown in Fig. ^ for J t — 1. The 
agreement with the above estimate for the 5 = 1/2 step is excellent. 

Another feature that can be fairly easily understood is the relative heights of the tips 
of the peaks in Fig. [| for states with 5 = 1/q. These states consist of one vortex every 
q plaquettes in the ground state configuration. The lowest energy domain wall consists of 
having one spacing that is q — 1 rather than q between two consecutive vortices. This should 



cost an interaction energy of about r ~ exp(— q/lo) with l £m of Eq. ( |TTD and indeed the 
tips of the peaks in Fig. for states with 5 = l/q fit this relationship quite well. 

The low temperature transverse resistance should be proportional to the thermal activa- 
tion rate of domain walls: Rt ~ exp(— 2Jefe//csT) (with e& replaced with e v for the Meissner 
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state). This facilitates direct comparison with experiment |J. 

IV. INCOMMENSURATE STATES AND PINNING-DEPINNING TRANSITION 

We now discuss the incommensurate states and superconducting-normal transitions in 
the transverse direction in these states. An incommensurate state is a state for which S is 
an irrational number. For J t = 0, the ground state has 7, = 7^ = and the aj are just 
pseudo- or slave variables determined by the constraint Eq.(§): 

oij = 27r/j + a - 2tt Y^~! n i- ( 12 ) 

The average vortex density (rij) is /; screening currents are absent, cto in Eq. (|T2|) is 
arbitrary; the a's are "unpinned" for all /. The system is simply two uncoupled ID XY 
chains (or rigid rotor model), so that the correlation length £ = 1/ksT (see for instance 
Ref. [Q). The system is superconducting at zero temperature along the ladder, but not 
in the transverse direction. As Jt rises above zero we observe a distinct difference between 
the system at rational and irrational values of /. For / rational, the a's become pinned 
for J t > (p(a) is a finite sum of delta functions) and the ladder is superconducting in 
both the longitudinal and transverse directions at zero temperature. For the / = case, the 
correlation length drops from the rigid rotor value of £ = 1 /k^T to the value found above in 
Eq.(|lTD like T- 1 ' 2 . 

The behavior for irrational / is illustrated in the following for the state with 5 = a g , 
where a g = (3 — v / 5)/2 (one minus the Golden Mean). Fig. |] displays p(a) for several 
different J t at H = a g . We see that the zero- frequency phonon mode (the smoothness of 
p(d)) persists for small Jt > until a critical value J t c (/) ~ 0.5 where the a's become 
pinned and the ladder becomes superconducting in the transverse direction. In the DSG, 
the pinning transition of this state coincides with the devil's staircase of Fig. |3| becoming 
complete |||| (If the a/s are pinned in this state, then all incommensurate states should 
be pinned). ^From Fig. |B| we see that this transition is a transition where the integration 
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measure p(a) in the partition function goes from a set of measure one on [— tt, tt) to a set 
of measure zero (at T — > 0). As one approaches the transition from J t < J£, p(a) — > 
continuously almost everywhere in [— tt, tt). As an order parameter for the transition, we use 
p(— 7r). This is shown in Figure ^ at a number of temperatures approaching T = 0. At the 
same J t the correlation length drops from the rigid rotor value of £ = l/lkgT) continuously 
towards a value of a few tens of lattice constants. This is shown in Figure |8|. For J t < J£ 
£ = l/(fcsT) is limited only by the temperature of the system. This is also true of the higher 
order correlation lengths, [In |Ao/A n ] _1 , n > 1 in the unpinned phase which also diverge like 
l/{k B T). 

Following standard scaling arguments one would expect p(— tt) to scale according to 

where j = J t c — Jt and £ ~ j~ u for T = while £ is a function of j and T for T > 0. For 
T — > 0, £ — > oo in the unpinned phase and p(— 7r) should be independent of £ implying 
/(y) — > constant as y —>■ oo so that p(— 7r) ~ j^. In the opposite extreme, with T fixed and 
j — > 0, p(— 7r) is dominated by the finite £ so that /(x) — » l/x' 3 as a; — * and p(— 7r) ~ ^^/ u . 
This allows us to define a more convenient scaling function /(a;) = l/x^f(x) so that 

p(-vr)(j,o = r /3/ v(je 1/l/ ). 

The results of this scaling, shown in the inset of Figure ^|, gives (3/v = 0.0926 ± 0.0009 and 
l/u = 0.352 ± 0.001, or /3 = 0.263 ± 0.003 and v = 2.841 ± 0.008. 
A similar scaling can be applied to £: 

aj,T) = T- 1 (a + b/\nT)g( J e / n 

where we have included the correction to scaling term b/ InT which causes a slight improve- 
ment to the fit away from J t c . (A power law correction was also considered, but did not have 
a statistically significant coefficient.) The resulting scaling collapse is shown in the inset of 
Fig. ^. This does not really give us any new information, but does act as a check on the 
information from the scaling of p(— 7r). 
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The pinning transition of the incommensurate states can be also studied using Eq. 
which are equivalent to the two-dimensional map: 

sin7j- + i = sin7j — J t sin a,, (13a) 
a j+1 = a,j - 27 j+1 + 2tt(/ - n j+1 ). (13b) 

Vortices are required in order to keep ctj in (— tt,tt]. Every trajectory of this map is a 
zero-temperature metastable state of the ladder. At Jt — the orbits of ( pT3|) for the 
incommensurate states will fill in a straight line 7 = in the 7a plane. As J t increases, 
these Kolmogorov-Arnold-Moser (KAM) orbits become deformed but remain smooth with 
the energy remaining independent of ao- Once Jt exceeds a critical value the KAM curve 
breaks down into a Cantor set of measure zero, as shown in Figure |9|. The disappearance of 
the zero-frequency phonon mode for irrational H's at finite small Jf(/) is equivalent to the 
breakdown of the Kolmogorov-Arnold-Moser (KAM) trajectories of the map ||15|| . 

However, we have been unable to find any connection between the exponents we found 
(describing the approach at H = a g to J t c and T = 0) to those found in ]n| for the approach 
to H = a g at J£ and T = for the standard map. If one expands the sinTj ~ jj in ( |13|) and 
redefine jj — ► 7^ — nf one obtains the standard map, studied by several authors including 
If one compares the critical J t c to the equivalent value for the standard map (J t = k c /2 



in the language of ]T5| ) one finds that they are identical. This is somewhat surprising. One 
might expect the two problems to be in the same universality class with related exponents 
but the location of the critical point normally depends on the details of the problem. 



V. EFFECT OF A CURRENT 

We now turn to the subject of critical currents along the ladder. One can obtain an 
estimate for the critical current by performing a perturbation expansion (i.e. {nj} remain 
fixed) around the ground state and imposing the current constraint of 

sinTj + sin 7^ = /. (14) 
11 



Letting 5jj and Satj be the change of jj and aij in the current carrying state, an expansion 
to first order gives 

cos7j_i5aj_i + cosjj5aj + i 
— (cos 7.,- + cos7j_i + 2 cosOj)5aj = 0, (15a) 

hj = - 5 % + ^")- ( 15b ) 

Eq. ( |15a| ) is in the form of G ■ 8a = 0. If detG 7^ 0, then 8ctj = and ^ = 1/2 cos 7^. 
The critical current can be estimated by the requirement that the 7, do not pass through 
7r/2, which gives I c = 2{n/2 — 7 max ) cos7 max , where 7 max = maxy^). In all ground states 
we examined for J t = 1.0, commensurate and incommensurate, we found that 7 max < vr/2, 
implying a finite critical current for all /. 

The presence of a current can also have an effect analogous to weakening the transverse 
coupling J t . One can get an approximate phase diagram in the I-f plane, valid for very 
low temperature and / < I c . To do this, rather than integrate out the f]i in Eq.(^|), one can 
substitute the current constraint Eq.(|TJ|) and get the "partition" function for the current 
carrying states as 

N r r / 1 

Z = \{ da^exp {(3J t cosai) cosh 2KJcos 2 [(a^i - a<)/2 + irf] - {I /J) 2 /A . (16) 

i J-TX L J 

Again, this is in the form of a product of transfer matrices and can be easily solved by the 
same techniques as before (see the appendix). The resulting phase diagram is shown in 
Figure [10] (compare to Fig. 

Comparing this diagram to Fig. ^ one can see that the increasing / in Fig. [H] has the same 
effect as decreasing J t in Fig. |j. This opens the possibility of experimentally studying things 
like the pinning-depinning transition experimentally. In this case, the pinning-depinning 
transition would be driven by the longitudinal current, and a measurement in the transverse 
direction would be used to probe the coherence of the two chains. 

This phase diagram in Fig. |TD| should be taken too literally however as once the critical 
current for the system is exceeded, the current constraint of Eq. (0) is no longer valid and 
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the system can switch between metastable states or may change continuously between a 
very large number of states. 



VI. CONCLUSION AND ACKNOWLEDGMENTS 

In conclusion, we have studied the equilibrium behavior of a Josephson junction ladder 
in a magnetic field in the absence of charging effects. Screening currents play an important 
role in this system, resulting in the spatial periodicity of the ground state climbing a devil's 
staircase as a function of /. Incommensurate states undergo a superconducting-normal tran- 
sition in the transverse direction as J t is increased, so that for J t > J£ « 0.5 the ladder is 
superconducting in both the longitudinal and transverse directions for all /. The critical 
current along the ladder is found to be finite for all /. Finally, although in one dimension 
there is no phase transition and long range order at finite temperature, our study showed 
that the correlation lengths in vortex states are extremely long for reasonably low tempera- 
tures. Thus one could test experimentally the predictions for the vortex configuration by, for 
instance, direct imaging via a scanning Hall-probe microscope or measuring the fractional 
giant Shapiro steps ]IB| . 

We thank Sue Coppersmith, Xinsheng Ling, and Qian Niu for valuable discussions. 

APPENDIX A: NUMERICAL SOLUTION TO INTEGRAL EIGENVALUE 

EQUATIONS 

To find the eigenvalues of the transfer matrix (0), we used the Nystrom method with 
n-point Gauss-Legendre quadrature ||17|| . The resulting matrix eigenvalue problem was then 



solved using the routines in LAPACK. It is worth describing this method here as it is fairly 
simple. We should also note that this method is quite fast, compared to methods such as 
the effective potential methods used by Mazo et al. in Ref. H. Indeed one can generate 



phase diagrams such as those in Fig. |3| in a few minutes or in Fig. |TD| in a hour or so of 
computer time (on an SGI workstation). 
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The transfer matrix P is not symmetric (except when / = or 1/2), and we will therefore 
require both right ip n and left ip n eigenvectors for the calculation of correlation functions. 
The eigenvalues and left and right eigenvectors are defined by the relations 

da'P(a,a')ijj n (a') = \ n ip n (a) (|A n | > |A n+ i|) 

-7T 

daip n (a)P(a,a') = \ n ip n (a'). (Al) 

The left and right eigenvectors are orthonormal (/ da ip n (a)ip n i (a) = 5 nn r) and complete (In 
all cases examined, the eigenvalues, though possibly equal in magnitude as in a complex 
conjugate pair, were nondegenerate). Note that the elements of the transfer matrix P are 
real and positive, so that the largest eigenvalue is real, positive and nondegenerate. Also, 
any complex eigenvalues will occur in complex conjugate pairs so that the partition function 
remains real and positive for any N. 

The Nystrom method requires the choice of some approximate quadrature rule: 

/ n n 
da'P(a,a')ip n (a') ~ ^2wjP(a, a'^ n a'^ (A2) 
j i 

Here the set {wj} are the weights of the quadrature rule, while the n points {a'j} are the 
abscissas. Since the solution method involves 0(n 3 ) operations it is worthwhile to use a 
high-order quadrature rule. For smooth, nonsingular problems Gaussian quadrature tends 
to be the best. Elements of the transfer matrix become increasingly singular as T — > 0, which 
produces computational difficulties. (For instance, a free energy per plaquette of 2.836 at 
kbT/ J = 0.004 corresponds to Ao ~ 10 308 , or close to overflow on the Silicon Graphics Indigo 
workstation used). Despite this problem, we found n-point Gauss-Legendre quadrature (n 
was typically 150 to 500) to give quite reasonable results. The largest errors tend to occur 
in the location of the edge of the steps of the devil's staircase, and even these errors are far 
too small to be visible on the plots. The lowest temperature achieved was k b T '/ J = 0.0007, 
which we believe is sufficient to characterize the states of the zero temperature system (to 
achieve this temperature, a constant was added to the energy so that the ground state energy 
was near zero, thus helping to correct for the overflow problem). 
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Equation is a standard eigenvalue equation 

P • = X^j (A3) 

which was solved using the routines in LAPACK (LAPACK library routines are provided 
on most workstation based systems and source code can be obtained free from netlib [0). 
When making use of the eigenvectors (such as in the calculation of correlation functions) 
at points not included in the quadrature points one should make use of Eq. ( |A2|) in order to 
maintain the full accuracy of the solution. In addition, in order to keep track of the weights, 
it is useful to symmetrize the weighting: defining the diagonal matrix D 1 / 2 = diag(yjwj). 
then the eigenvalue equation becomes 

( D l/2p D l/2) . (D l/2^ = A(D 1 /V), 

which is in the form of a symmetric eigenvalue problem for / = 0, 1/2. 
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FIGURES 
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FIG. 1. The Josephson junction ladder is formed by the arrangement of the superconducting 
islands. The field H is out of the page and the arrows indicate the direction of the gauge invariant 
phase differences. 
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FIG. 2. Correlation function for a three periodic state (H = 1/3). The three branches cor- 
respond to (e'H-asn)^ ( e i(a -a 3n+1 )^ and ^i(ao-«3n+ 2 )^ with n = o,l,2,---. The dotted line 

indicates the value of (e ta i) 2 . 
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FIG. 3. S = Arg(X 1 )/27r versus / for k B T/J = 0.002 and (a) J t = 0.4, (b) J t = 0.6, and (c) 
Jt = 1.0. 
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FIG. 4. Periodicity S phase diagram for low order rationals. At finite temperature, the step 
edges are slightly rounded. The phase boundaries for a step at E = p/q enclose all E satisfying 
|E —p/q\ < e, with dashed lines indicating an e of 10~ 6 and solid lines indicating an e of 0.002 at 
k B T/J = 0.0007. 
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FIG. 5. Effective Ising coupling as a function of / for J t = 1. Inset: Statistical error for 2e& in 
the fit versus /. 




FIG. 6. p(a) = tp^(a)ip^(a) versus a and S = 0.381 966 011 ••• , and for (a) J t = 0.4, (b) 
J t = 0.486, (c) J t = 0.59 at k B T = 0.0007J and (d) J t = 0.7 at k B T = 0.001 J. Note the smaller 
scale for the upper plots. 
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FIG. 7. p(-7r) versus J t for k B T = 0.012J, 0.008J, 0.004J, 0.002J, 0.001J, and 0.0007J. Also 
shown is the extrapolation to T = 0, 0.192257(J t c — Jt) 13 ■ Inset: Scaling collapse of same data 
shown in main plot. 
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FIG. 8. Correlation length £ versus J t for /c B T = 0.012J, 0.008J, 0.004J, 0.002J, 0.001J, 
and 0.0007J. Inset: Scaling collapse of same data shown in main plot (a = 0.793 ± 0.006 and 
b = -0.39 ±0.04). 
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FIG. 9. KAM curves constructed from plotting each (aj,7j) as a point in the plane. The state 
shown is one with periodicity 28657/75025 w 0.381966 01 at / = 28657/75025. Below the critical 
coupling J t c , the points fill in a smooth curve (top panels) whereas above J t c the points form a zero 
measure Cantor set. Insets: Histograms of the {<x,}, which can be compared to the phase density 
p(a) calculated from the transfer matrix 
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FIG. 10. Periodicity as a function of current / and filling factor / for J t = 1.0. At finite 
temperature, the step edges are slightly rounded. The phase boundaries for a step at S = p/q 
enclose all S satisfying |E —p/q\ < e, with dashed lines indicating an e of 10~ 6 and solid lines 
indicating an e of 0.002 at ksT/J = 0.004. Note that some of the "tongues", such as the one 
for S = 1/4 come to a constriction at what appears to be the critical current for that state. The 
section above the constriction is, therefore, above the critical current and should not be taken too 
seriously. 
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